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Monte Carlo simulations are performed to determine the critical percolation threshold for inter- 
penetrating square objects in two dimensions and cubic objects in three dimensions. Simulations are 
performed for two cases: (i) objects whose edges are aligned parallel to one another and (ii) randomly 
oriented objects. For squares whose edges are aligned, the critical area fraction at the percolation 
threshold cj> c = 0.6666 ± 0.0004, while for randomly oriented squares (f> c = 0.6254 ± 0.0002, 6% 
smaller. For cubes whose edges are aligned, the critical volume fraction at the percolation threshold 
4> c = 0.2773 ± 0.0002, while for randomly oriented cubes C = 0.2236 ± 0.0002, 24% smaller. 



I. INTRODUCTION 

Lattice percolation is often used for the statistical mod- 
eling of transport in porous media Jl]-^| . The requirement 
that sites, and the bonds between them, be constrained 
to a fixed lattice may not, however, be an appropriate 
model for natural porous media Jl]-|| . The characteristics 
of site and bond percolation potentially limit their appli- 
cability to modeling of natural phenomena such as oil 
or groundwater flow and extraction of melt from super- 
solidus regions deep inside a planetary body. 

Continuum percolation offers two advantages for de- 
scribing porous media: 

(i) The objects which form clusters are not restricted 
to points on a fixed lattice; they can be placed 
anywhere within the volume studied and either 
be barred from interpenetration, or allowed to in- 
terpenetrate, i.e., they can have either "hard" or 
"soft" cores B. Because of the freedom of place- 
ment inside the system, the connections between 
soft core objects can range from very small to very 
large, depending upon the extent of interpenetra- 
tion. 

(ii) The objects can be of any shape. In two dimensions 
the continuum percolation of discs is often investi- 
gated |^,||- In most studies of continuum perco- 
lation in three dimensions, spheres are used as the 
objects, leading to the "swiss-cheese" nomenclature 
for continuum percolation Other frequently- 
used shapes are rods and ellipsoids of revolution || . 
In a few cases the continuum percolation of cubes 
has been considered ||. 

Here we determine the threshold for continuum perco- 
lation of soft core squares in two dimensions (2-D) and 
cubes in three dimensions (3-D) whose edges are aligned 
parallel to, or oriented at random angles to, the axes of 
the system. Continuum percolation is believed to belong 
to the same universality class as site and bond percolation 



j7|,||; once we have determined the continuum threshold 
for an object of a specific shape we can apply many of the 
characteristics of site and bond percolation, e.g. critical 
exponents, to describe the continuum percolation cluster. 



II. METHODS 

We construct 2-D and 3-D Monte Carlo simulations 
for the determination of the percolation threshold based 
upon the Leath method || and the methods Lorenz and 
Ziff |n]] used in their study of the continuum percolation 
of spheres. 

We perform 2-D simulations with squares whose edges 
are of unit length. In 3-D we perform simulations with 
cubes of three different edge lengths: l/VS, 0.75, and 1.0; 
the length L of each axis of the simulation box in the 2-D 
simulations is 301, and in 3-D is 101 in simulations using 
cubes with edge lengths of l/\/3 and 0.75, and is 161 in 
simulations with the unit cube Q. We subdivide the 
system into either a 2-D or 3-D grid of unit area squares 
or unit volume cubes; an illustrative 12 x 12 2-D version 
of our system is shown in Fig. [I]. 

The cluster begins in the center grid volume and ob- 
jects are added to it based upon a Poisson distribution 
centered about the average number of objects per unit 
area or volume, N/L d , chosen for the simulation of di- 
mension d JlOf . The product of this value and the indi- 
vidual object's area or volume, v, is the reduced number 
density 
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(i) 



If the number of objects, n, generated from the Poisson 
distribution is nonzero, then these n objects are placed 
at random locations inside the grid volume. To fix the 
orientation of each individual object, random numbers 
are generated to determine the one angle of rotation in 
the 2-D simulations and the three Euler angles of rota- 
tion in the 3-D simulations; these angles vary between 
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and 2tt. The locations of the 4 (2-D) or 8 (3-D) cor- 
ners and the center of each object are stored in a data 
structure, along with a flag indicating that the grid area 
(2-D) or volume (3-D) was visited and populated during 
the realization. The nearest-neighbors and next-nearest- 
neighbors of this grid area or volume are then populated 
in a similar manner and the intersection between squares 
or cubes is tested. 

To determine whether two squares intersect we choose 
one square in the cluster as the reference square and an- 
other square in the simulation as the test square. We use 
an algorithm for the intersection of two lines |L2| to test 
if any of the 4 edges of the reference square interesect 
the 4 edges of the test square. To determine if two cubes 
intersect we choose one cube in the cluster as the refer- 
ence cube and another in the simulation as the test cube. 
We use an algorithm for the intersection of a line and a 
facet jL3| to test if any of the edges of the reference cube 
intersect the faces of the test cube. In this algorithm, 
the location of each of the 12 edges of the reference cube 
are compared to the location of the 12 triangular facets 
that describe the locations of all faces on the test cube 
using the corners and diagonal of each face. If the test 
object intersects the reference object, it is added to the 
growing cluster. This process is repeated for each new 
square or cube added to the system until the cluster can 
no longer grow. Intersections between squares or cubes 
in grid areas or volumes up to two units away can oc- 
cur for edge lengths 0.75 and 1, as exemplified in 2-D 
for squares of unit edge (Fig. [l]) , but cubes of edge 1 / \/3 
can only interesect if they are in the same or neighbor- 
ing volumes of the grid, which reduces the number of 
grid volumes that must be checked for cube overlap in 
simulations with cubes of this smallest size. 

The cumulative distribution of cluster sizes is calcu- 
lated from the cluster size of each realization, s, by bin- 
ning the cluster sizes such that all bins in the range of 
2° to 2 S+1 — 1 are incremented by 1. In order to esti- 
mate finite-size effects of the simulation, objects in each 
cluster are tested to determine whether they touch the 
edge of the simulation. If so, the cluster size is compared 
against the smallest cluster size in previous realizations 
that touch the edge and the smaller value is stored. Bins 
of size greater than the smallest cluster that touched the 
edge of the system are not used in the determination of 
the percolation threshold. 

At the end of the simulation, the value of each bin 
is divided by the number of realizations, from 10,000 to 
50,000, to yield the probability of achieving a cluster of 
size s, P(s\rj), for a given value of r\. Power law behavior 
of the probability as a function of the bin size is inter- 
preted to indicate the critical percolation theshold r; = r/ c 
[^,1). To accurately determine the theshold, we follow the 
techniques of Ref. |ici|j . The probability of generating a 
cluster size s at a specified r\ is 

P(s\t,)~ As 2 -f[(rj- Vc )s^}, (2) 

where both r and er are universal exponents and A is 



a non- universal constant. In 2-D the values of these 
exponents are 187/91 and 36/91 [Q, respectively. In 
3-D the values of r and a are 2.18906 ± 0.00006 and 
0.4522 ± 0.0008, respectively [||). Near the percolation 
threshold the scaling function fix) can be expanded in a 
Taylor series 

fix) = 1 + Bx + 0(x 2 ). (3) 

Combining Eqs. (||) and (^), 

P(s\ri) S T - 2 ^ A + AB{ V -T lc ) S a + ■■■ (4) 

which demonstrates that P(s\r/)s T ~ 2 becomes constant 
at the percolation threshold as s becomes asymptotically 
large. 



III. RESULTS 

The percolation threshold can be expressed either as 
the critical reduced number density, rj c , or the critical 
area (or volume) fraction, <fi c , which are related to each 
other by |p| 

4> c = \-e~^. (5) 

A. Two Dimensions 

For square objects aligned parallel to each other in the 
2-D system we find 

r) c = 1.098 ±0.001 (6a) 

or, from (||), 

c/) c = 0.6666 ± 0.0004 (6b) 

(see Fig. ||a). Our value of <\> c is within the error bars of 
two previous determinations by Monte Carlo techniques, 
where C = 0.668 ± 0.003 JT§] and <j) c = 0.65 ± 0.02 
Jl7| . However our determination of <f> c is slightly lower 
than that calulated by || whose Monte Carlo simula- 
tions produced <fi c = 0.6753 ± 0.0008, and whose ap- 
plication of the direct-connectedness expansion method 
yielded (j) c — 0.6912. In contrast, our value of 4> c is sig- 
nificantly higher than the experimental one of |18) whose 
average for 9 trials is <fi c — 0.613 ± 0.013. 

We find that for randomly oriented square objects in 
2-D 

? ?c = 0.9819 ± 0.0006 (7a) 

or 

C = 0.6254 ±0.0002 (7b) 
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(see Fig. |2|b). This is the first determination of these 
values for randomly oriented squares. 

These values for the continuum percolation threshold 
for aligned and randomly oriented squares are lower than 
for discs, <f> c — 0.67 [^9|, by a maximum of sa 0.5% for 
aligned squares and w 7% for randomly oriented squares. 
We attribute the significant difference in </> c between discs 
and randomly oriented squares to the possibility of ran- 
domly oriented squares intersecting other squares whose 
centers are located at distances up to the diagonal length 
of the square (see Fig. [I]), whereas two discs can only in- 
tersect if their centers are no further than one diameter 
away from each other. The similarity of <j) c for aligned 
squares and discs may occur because both objects can 
only possibly intersect other objects whose centers are 
separated at most by either the edge length of the square 
or the diameter of the disc. 



B. Three Dimensions 

For cubic objects aligned parallel to each other in the 
3-D system, 

r/ c = 0.3248 ± 0.0003 (8a) 

or 

<j) c = 0.2773 ±0.0002 (8b) 

(see Fig. ||a). The precision of this result is greater than 
the most precise previous determination <f> c — 0.280 ± 
0.005 ||. The critical volume fraction is significantly less 
when cubic objects are allowed to have random orienta- 
tions, 

r, c = 0.2531 ± 0.0003 (9a) 

or 

<j) c = 0.2236 ± 0.0002 (9b) 

(see Fig. ||b). The result for randomly oriented cubes 
is the same for cubes of edge-length l/\/3, 0.75, and 1. 
Thus, as expected, the percolation threshold is indepen- 
dent of the cube and system size used. This value is the 
first determination of <f> c for the continuum percolation 
of randomly oriented cubes. 

Comparision of the critical volume at the percolation 
threshold for aligned cubes with that determined for 
spheres, cj) c = 0.289573 ± 0.000002 @, demonstrates 
that the difference in shape between spheres and cubes 
affects 4> c by « 4% || . Allowing cubes to randomly orient 
lowers 4> c ~ 23%. The difference between the randomly 
oriented cubes and spheres is due to the same process as 
discussed above for discs and squares, but in this case it 
is the greater length of the body diagonals of cubes com- 
pared to the diameter of spheres or the edge length of 
aligned cubes that enhances the probability of connect- 
edness for randomly oriented cubes at any given volume 
fraction. 



IV. DISCUSSION 

The continuum percolation threshold can be predicted 
with excluded volume theory |^, |20| , ^2]| : 

N c V ex = B c , (fO) 

where N c is the critical density of objects [p3[ , V cx is their 
excluded area or volume, and B c is the average number of 
bonds per object |2l}| . Originally, B c was thought to be 
one constant for all parallel (i.e., not randomly oriented) 
convex objects in 2-D and another constant in 3-D 
but later B c was determined to be different for spheres 
and for cubes in 3-D ||. The excluded area for discs 
and aligned squares of unit area is 4 and for randomly 
oriented unit squares 4.084 pc| . For both spheres and 
aligned cubes V ex is equal to 8 times their volume in 3-D. 
For randomly oriented cubes Vex is 11 times their vol- 
ume |24|]. Calculated values of B c in 2-D and 3-D are 
presented in Table I. 

In 2-D we determine B c = 4.39 ± 0.01 for aligned 
squares and for randomly oriented squares B c = 4.01 ± 
0.01. The value for aligned squares is similar to that 
originally proposed by Balberg for discs and squares, 
B c — 4.5 ±0.1 |||] and to the values for discs and squares 
calculated from Monte Carlo simulations: B c = 4.43 (H]] 
and B c = 4.5 ±0.1 ||, respectively. On the other hand 
our value for B c is somewhat lower than calculated by a 
series expansion technique, B c — 4.7 |3|. Thus we confirm 
that B c has the same value, within error, for discs and 
for aligned squares in 2-D. However, a different behavior 
is observed when squares are randomly oriented and B c 
drops by ~ 10%, consistent with expectations that non- 
parallel objects should exhibit a lower B c than parallel 
objects 

Our Monte Carlo simulations yield B c = 2.59 ± 0.01 
for aligned cubes, as is expected because of the agree- 
ment between our estimate of the percolation threshold 
and previous estimates. For randomly oriented cubes 
B c = 2.78 ± 0.01, within error of B c for the continuum 
percolation of spheres This result is surprising in 
light of the observation that B c for randomly oriented 
squares is » 10% below aligned squares and is contrary 
to expectations that randomly oriented objects should 
have lower values of B c than aligned ones p3[ . However, 
B c for randomly oriented cubes does not exceed the lim- 
iting value predicted by the excluded volume theory of 
the continuum percolation threshold pl| . 

Our results confirm previous research demonstrating 
the effect of object shape on the threshold for contin- 
uum percolation. We furthermore find that the incor- 
poration of random orientations of objects in continuum 
percolation simulations significantly affects the percola- 
tion threshold. Most of these effects are predicted by the 
application of excluded volume theory to the calculation 
of the percolation threshold but, surprisingly, allowing 
squares to randomly orient decreases B c in 2-D whereas 
random orientation of cubes increases B c in 3-D. 
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FIG. 1. Two-dimensional, 12 x 12 example of a percolation cluster of unit- length square objects (thick lines) for the case of 
randomly oriented squares. The system is divided by a series of grid lines (thin lines) that create unit areas in this 2-D system. 
Note that the upper two objects of the cluster in the center of the system intersect each other even though their centers are 
placed in next-nearest-neighbor areas of the grid. The real 2-D and 3-D systems of our study are much larger than this system. 
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(b) Randomly Oriented Squares 
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FIG. 2. Power-law scaled plots for determination of percolation threshold for squares of unit size in a 301 x 301 system based 
upon 50,000 realizations at each area fraction. At the threshold rj — rj c , P(s\rj)s T ~ 2 is independent of s" , which allows for 
accurate determination of rj c , which is related to cf> c by Eq. (a) Squares whose edges are aligned parallel to each other, for 
which case we estimate r\ c — 1.098 ± 0.001, so (f> c = 0.6666 ± 0.0004 by Eq. [E|. (b) Squares that are randomly oriented as shown 
in Fig. 1, for which we estimate r) c = 0.9819 ± 0.0006, so </> c = 0.6254 ± 0.0002. 



7 



1.00 $ 



(a) Aligned Cubes 




A Volume Fraction = 0.2776 
O Volume Fraction = 0.2774 
V Volume Fraction = 0.2772 



20 



30 



40 



50 



1.25 



1.20 - 



1.05 -V 



1.00 $ 



(b) Randomly Oriented Cubes 
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FIG. 3. Power-law scaled plots for determination of percolation threshold for cubes of unit size in a 161 x 161 x 161 3-D 
system based upon 50,000 realizations for each volume fraction. At the threshold r/ — r) c , P(s\rj)s T ~ 2 is independent of a" , 
which allows for accurate determination of rj c , which is related to <f> c by Eq. | (a) Cube; whose faces are aligned parallel to 
each other, for which case we estimate r\ c = 0.3248 ± 0.0003, so by % <t> c = 0.2773 ± 0.0002. (b) Cubes that are randomly 
oriented, for which we estimate r) c = 0.2531 ± 0.0003, so C = 0.2236 ± 0.0002. 
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